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Starting with stochastic rate equations for the fundamental interactions between microbes and 
their viruses, we derive a mean-field theory for the population dynamics of microbe-virus systems, 
including the effects of lysogeny. In the absence of lysogeny, our model is a generalization of that 
proposed phenomenologically by Weitz and Dushoff. In the presence of lysogeny, we analyze the 
possible states of the system, identifying a novel limit cycle, which we interpret physically. To test 
the robustness of our mean field calculations to demographic fluctuations, we have compared our 
results with stochastic simulations using the Gillespie algorithm. Finally, we estimate the range of 
parameters that delineate the various steady states of our model. 
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I. INTRODUCTION 

Microbes and their viruses are the most genetically di- 
verse, abundant and widely distributed organisms across 
the planet [ll-Ql • Microbes are major contributors to the 
global biogeochemical cycles and catalyze the reactions 
that have over evolutionary time brought the Earth's sur- 
face to its present redox stateQ. Similarly, viruses, es- 
pecially in the oceans, manipulate marine communities 
through predation and horizontal gene transfer 0|, re- 
cycle nutrients and thus drive the biological pump which 
leads inter alia to the sequestration of carbon in the deep 
ocean [i-fTrl. 

It is being increasingly realized that the classical view 
of microbial viruses purely as predators is too limited. 
Many microbe-virus interactions are lysogenic, not lytic: 
upon infection, the viral genetic material is incorporated 
into the chromosome of the host, replicates with the 
host, and can be subsequently released, typically trig- 
gered by the stress response of the host to environmental 
change [ISj. As a result, viruses can transfer genes to 
and from bacteria, as well as being predators of them, 
so that the virosphere shoul d prop erly be recognized as 
a massive gene reservoir [TtI [l9l - l2]| . Thus there is a co- 
evolution of both communities, the effects of which are 
complex and far-reaching [l^, [13, l20l - [25| . even including 
the manipulation of bacterial mutation rates [1^. This 
nontrivial interaction between microbes and viruses has 
not gone unnoticed, with wide interest among biologists, 
ecologists and geologists dHUHlMl • 

These findings highlight the importance of consider- 
ing ecosystem dynamics within an evolutionary context. 
Conversely, evolution needs to be properly understood 
as arising from a spatially- resolved ecological context, as 
was first recognized by Wallace over 150 years ago|35j. 
That speciation, and adaptation in general, arises at 
a particular point in time and space has a number of 
deep consequences that have not yet been incorporated 
into current theory. First, it means that evolutionary 



dynamics proceeds by the propagation of fronts, result- 
ing in a complex and dynamical pattern of speciation, 
adaptation and genome divergence that reflects its in- 
trinsic dynamics and that of the heterogeneous and dy- 
namical environment 36 39] . Second, as fronts expand, 
there are only a few pioneer organisms at the leading 
edge, and so demographic fluctuations are much larger 
than in the bulk. Such fluctuations profoundly influ- 
ence the spatial structure of the populations, and during 
the last few years have been recognized to play a ma- 
jor role in population cyclesjiO] and even spatial pattern 
formation[4l|. Third, the existence of horizontal gene 
transfer and genome rearrangement processes is strongly 
coupled to spatial distribution. For example, it is known 
that the probability of conjugation events is dependent 
on the local density, being essentially one/per genera- 
tion in closely-packed biofilms, but an order of magni- 
tude smaller in planktonic culture [l^. Moreover, the 
mechanism of horizontal gene transfer is also dependent 
on the density, with viral-mediated transduction being 
the most relevant mechanism at low density. How these 
patterns of evolutionary dynamics and species distribu- 
tion play out is essentially unexplored. However, there 
have recently been the first reports of observations of the 
coupling between evolutionary and ecological timescales. 
In one such system (a predator-prey system realized in 
rotifer-algae interactions) , it has been demonstrated that 
rapid evolutionary dynamics is responsible for the un- 
usual phase-lag characteristics of the observed population 
oscillations [1^. Thus, rapid evolution is not only a ma- 
jor force for adaptation, but can have marked ecological 
consequences too. 

The complex interplay between evolution and the en- 
vironment is nowhere more important than in early life, 
where the key questions concern how life emerged from 
abiotic geochemistry. Early life experienced demanding 
environments, whose closest modern day correspondence 
might be deep ocean hydrothermal vents or hot springs. 
It is known that there are high occurrence of lysogens in 
both environments [1, 113], suggesting that microbe-phage 
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interactions might also be important in the early stages of 
life, with lysogens playing an important role as a reservoir 
of genes and perhaps even aiding in the stabilitization of 
early life populations through the limit cycle mechanism 
discussed in this paper. 

Our goal in this paper is to lay a theoretical foundation 
for describing the interplay between ecology and evolu- 
tion in the context of microbe-virus systems, as these 
are arguably amongst the most important and proba- 
bly the simplest of the complex systems in biology. The 
questions that will ultimately interest us are the evolu- 
tionary pressures that tune genetic switches governing 
the lysis-lysogeny decision, as well as the factors that 
shape pr ophage induction in response to environmen- 
tal stress |45l-t47|. Such a foundation must begin with a 
proper account of the population dynamics itself, before 
coupling in detail to other levels of description involving 
genome dynamics, for example. Thus, we have chosen to 
focus in the present paper on the dynamics of microbe- 
virus systems, taking full account of both of the major 
viral pathways. In this paper, we are not specific about 
whether we are dealing with bacterial or archaeal viruses, 
but because most of the experiments to date are carried 
out on bacteria, we have tended to identify the microbes 
as bacteria and the viruses as phages, even though this 
is not required by the mathematics. 

We are now ready to introduce the specific problem 
that we treat in this paper. Upon phage infection, there 
are two pathways awaiting the host bacterium^]. In the 
first pathway — lysis — the bacteriophage produces a large 
number of copies of itself utilizing the bacterium's ge- 
netic material and molecular machinery. As a result, the 
bacterium ceases its metabolic function, and ruptures, 
releasing the newly-assembled bacteriophage inside. The 
other pathway is lysogeny. In this process, the intruder 
integrates its own DNA into the genome of the bacterium, 
enters a dormant stage and becomes a prophage. The in- 
fected bacterium is known as a lysogen — a relatively sta- 
ble state [49} ■ immune to superinfection from the same or 
related strains, and proceeding under normal replication 
life-cycles. The DNA of the bacteriophage is duplicated, 
along with that of the host during cell replication. The 
lysogenic state can be terminated by environmental stress 
such as starvation, pollution or ultraviolet irradiation, re- 
sulting in the process known as prophage induction: the 
exit of the prophage from the host genome, and sub- 
sequent lysis of the original bacterium and its bacterial 
descendants. 

We now discuss briefiy existing treatments of popula- 
tion dynamics in the context of microbe-virus systems. 
In 1977, Levin et al [s^l extended the celebrated Lotka- 
Volterra equations to model the dynamics between viru- 
lent phages and their victims, where only virulent phages 
are considered. A number of extensions have been pro- 
posed, extending the level of biological realism to in- 
clude such features as the time delay arising between 
infection and lysi s as well as the evolution of kinetic 
parameters fslmSj . In 2007, Weitz and DushofffSGj pro- 



posed another way to improve the classic predator-prey 
model. Their attempt was mainly based on the experi- 
mental observation that the ability of a bacteriophage to 
lyse hosts degrades when the bacteria approach their car- 
rying capacity (57, .59]. Adding a new term to account for 
the saturation of the infection of the bacteriophages, they 
obtained an interesting phase diagram in which the fate 
of the bacteria-phage community can depend on the ini- 
tial conditions. However, the new term is put in by hand, 
based on intuition which needs detailed mathematical 
support. Furthermore, they focused on virulent phages 
and excluded the temperate ones that elicit lysogeny, now 
regarded as essential to the survival of microbial commu- 
nities through fluctuating environments [isi [29l |33| . 

These works are based upon an ensemble-level de- 
scription of the community, as in the classic work on 
predator-prey svstems(60j. However, as is well-known 
[601, the simplest of these models fails to capture the 
intrinsic cyclical behavior of predator-prey populations 
despite apparently incorporating fully the basic inter- 
actions that should give rise to cycles. This paradox 
was resolved by the important work of McKane and 
Newman [i^l, who showed that cyclical effects could only 
be captured at the level of an individual-level model, and 
arose through the amplification of demographic noise. 
Their work showed how the conventional ensemble-level 
equations for predator-prey systems arose as the mean 
field limit of the appropriate statistical field theory, with 
the essential effects of demographic noise entering the 
analysis as one-loop corrections to mean field theory, in 
an inverse population size expansion. These effects can 
also be treated in a slightly more convenient formalism 
using path integrals jGlj . The literature also does not have 
an explicit representation of lysogeny as it modifies the 
population dynamics of both host and phage. 

The purpose of this paper is to provide a detailed the- 
ory of the population dynamics for host-phage commu- 
nities. In contrast to earlier work, we pose the problem 
at the microscopic level, working with an individual-level 
model of bacteria and phage. From this fundamental de- 
scription, we are able to derive the usual community level 
description analogous to Lotka-Volterra equations from a 
mean field theory. Our results encompass both virulent 
phages, such as those in Weitz and Dushoff's work[56j, 
and lysogenic phages which have not been studied math- 
ematically up to now. 

This paper is organized as follows. In Section HIl as a 
preliminary exercise, to present the technique, we treat 
a lysis-only model, in which we derive a set of dynami- 
cal equations roughly in the same form as in Weitz and 
Dushoff's paper[56| except for an additional parameter, 
which generally results in a relatively unimportant shift 
in the phase diagram. In the full lysogeny-lysis model, 
presented in Section Hill we develop the formalism for the 
community of hosts and phages, including both lysis and 
lysogeny. Interestingly, we find that for certain combina- 
tion of parameters, the community exhibits a limit cycle 
for all the species in the phase space, even at the level of 
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mean field theory. In order to interpret the correspond- 
ing range of parameters in a useful way for experimental 
observations, we map the parameters to rates in chemical 
reactions. In order to explore the robustness of our re- 
sults, we demonstrate in Section HVl that the correspond- 
ing limit cycle arises also in stochastic simulations with 
the Gillespie algorithm. Finally, in Section|V]we estimate 
the feasibility of verifying our predictions in laboratory 
experiments. 



for hosts and W units for phages. These units will be re- 
ferred to as the host layer and phage layer, respectively. 
In the host layer, each unit is either occupied by one host 
or unoccupied, i.e. an empty site E. The total number 
of empty sites E is K — in. We construct the phage layer 
in a similar manner and denote the empty sites there by 
(/) although the phage population is not confined actually. 

The population dynamics of the system can thus be 
modeled as arising from the following six microscopic 
events (Table IJ): 



II. LYSIS-ONLY MODEL 



TABLE I: Microscopic events in the lysis-only model. 



Derivation of the population dynamics from an 
individual-level model 



In this section, we adapt the classic predator-prey 
model to the host-phage communities from a microscopic 
or individual-level model. For simplicity, we first focus on 
two-component competition, where lysogens are excluded 
in spite of their biological importance. Hence, we are 
considering virulent phages and their hosts. Following 
the procedure given by McKane and Newmanjioj. we de- 
rive the population dynamics for the host-phage system, 
which Weitz and Dushoff,56j had written down based on 
intuition. Here we work at the level of mean field the- 
ory, and we do not, in this paper, include the extension 
necessary for representing spatial degrees of freedom. In 
our model, the host-phage dynamics differentiates itself 
from the classic predator-prey model in two ways: (1) 
only the host population is restricted by carrying capac- 
ity due to resource limitation and (2) the lysis of one 
host releases a particular number of phages (for exam- 
ple, about 100 replicates for lambda phage |48{). instead 
of only one predator in the classic predator-prey model. 
The above two points need to be accounted for carefully 
in the set up of the model, especially in the introduction 
and application of the carrying capacity, which will be 
explained explicitly as follows. 

In our host-phage community, we have only one species 
of host and one species of phage which preys upon the 
former. Let us label the hosts by A and phages by B, 
whose populations are m and n, respectively. The hosts, 
either heterotrophic or autotrophic, need to consume en- 
vironmental resources, which are renewable in every cy- 
cle, for survival and reproduction. All the environmental 
limitations on the hosts are abstracted into a maximal 
host population, which is denoted by the carrying capac- 
ity K. The phages, on the other hand, do not rely on the 
consumption of natural resources for maintenance once 
they are released into the environment. Thus, there is no 
such hard constraint on the phage population. Although 
phages are not restricted by K, we still introduce a vir- 
tual carrying capacity W for phages from dimensional 
considerations. It can be imagined that ^> oo so that 
no true carrying capacity is imposed on the phage pop- 
ulation. The carrying capacities can be better visualized 
if we conceive space to be equally divided into K units 



description 



symbols 



birth of host AE AA 

death of host due to longevity A E 

death of host due to crowding effect AA -4 AE 
host-phage interaction 

■ under good metabolism AEB A EEaB 

■ under poor metabolism 
death of phage 



AAB 4 EAI3B 



Here, 6, c, d, e, / and g are all constant rates. All the 
events above are written with constraints, with a nonlin- 
ear relation being incorporated automatically by adding 
empty sites E to the left of the arrows to reflect the re- 
striction of carrying capacity K. For example, the birth 
of the host is density dependent, which needs an empty 
site to accommodate the newly-born host. If no empty 
site is found, such an event can not happen. Since we con- 
sider only the mean field case, no spatial inhomogeneity 
is introduced. There is no concept of locality here, ei- 
ther. As long as an empty site is found, the newly-born 
host is permitted. The crowding effect describes the com- 
petition in survival for limited natural resources among 
hosts. No such crowding effect exists for phages, which is 
in line with our assumption that there is no true carrying 
capacity confining the phage population. The two events 
in host-phage interaction are carefully chosen to give a 
minimal model while encompassing reduced lysis when 
the host population is approaching its carrying capacity. 
In the events, a and (3 are the numbers of progeny for 
phage reproduction under good and poor metabolism, re- 
spectively. In biology, there are mainly two reasons which 
may account for reduced lysis. The first is the decrease 
in the phage reproduction number [STf, i.e. 



a > /3, 



(1) 



because phages need bacterial genetic materials, molecu- 
lar machinery and energy in the synthesis of their repli- 
cates. When the normal function of the host is down 
regulated, phage replication is correspondingly down 
shifted. Another reason is the reduced efficiency dur- 
ing phage infection, either in adsorption rates or viable 
infection, which leads to a diminishing of the infection 
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cvcle[57|. i.e. 



Thus, the time evolution for the population size is 



e>/. 



(2) 



It might seem as if the model is discrete in the repre- 
sentation of metabolism since we put in good and poor 
metabolism by hand. However, note that the actual 
metabolism of the community may be somewhere be- 
tween good and poor, i.e. a linearly combination, de- 
pending on the probability or fraction to enter either 
event. Finally, the death of the phage may be caused 
by protein cleavage. 

The time evolution of the whole community is accessed 
by random sampling. In each time step, we have a proba- 
bility /i to draw units in the host layer and a probability 
V to draw units in the phage layer. In the host layer, 
we may draw either one unit with probability w or two 
units with probability \ —uj. In the phage layer, only one 
unit is drawn. If a combination not listed in Table U is 
drawn, such as EEB, nothing happens. Thus all we need 
to consider are the above events. Using simple combina- 
torics, it is straightforward to obtain the probability for 
the combinations as follows (Table HI)): 



d (m) 
dt 



d{n) 
dt 



{T (m + 1, n|m, n)) — {T (m — 1, n\m, n)) 

— (T (to — l,n + a — 1|to,7i)) 

- (T(to - + - l|TO,n)) ; (4a) 

(a — 1) (T (to — l.n + a — l|TO,ri)) 
+ - 1) {T {m - l,n + /3 - l\m, n)) 



(T (to, n — 1|to, n)) . 



(4b) 



Here we have taken the mean field theory limit and 
neglected all the correlations and fluctuations. 

Omitting angle-brackets for simplicity, the equations 
for the evolution in population are 



dm 
'dt 

dn 



drr,m 



/ TO 

- (pmn ~ am — 



(5a) 
(5b) 



TABLE II: Probabilities for the combinations in the lysis-only 
model. 

combination probability 



where 



A 
AA 
AE 
AEB 

AAB 
B 



m(i - ^I'^Jl 



m (m — 1) 

K {K ~1) 
2m (K - m) 



fiiy (1 — ui) 



liv (1 — uj) 



K(K-l) 
2m {K — m) n 

K {K-1) W 
m (m — 1) n 

kIk-i)w 



w 



where the factor 2 accounts the equality in probability 
for events AE and EA, or AEB and EAB. 

Thus we obtain the transition matrices for each kind of 
variation in the population during each time step, such 
as (T (to -|- 1, n\m, n)), and further the evolution for the 
probability in the population with to hosts and n phages 
at time t P (to, n, t). The reader is referred to Appendix 
A for calculational details. 

The average of the population is given by summation 



(26 + d)M(l - v) (1 - w) 

r = : 

K 

2e^w (1 — Lo) 
^ " KW ' 
7 = a - 1; 

(cw + d{l - uj))fi (1 - ly) 



K 



W ' 
l_ 

2e' 

2ae' 



1 - 
1 - 



Considering Eq. ([T]) , we notice that Eq. ([6f 
the following relation 



< am < a„ < 1. 
Generally speaking, a™ ^ a„ unless 



(6a) 

(6b) 
(6c) 

(6d) 
(6e) 
(6f) 
(6g) 

yield 
(7) 



i) = mP (to, n, t) ; 

ran 

i) = nP (to, n, t) . 



(3a) 
(3b) 



(8) 



which implies that the reproduction numbers under good 
and poor metabolism are the same as in Weitz and 
DushofF's model. This concludes the derivation of the 
equations for population dynamics from the individual 
or microscopic level. 
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B. Results 



The last is the coexistence of hosts and phages 



In this section we explore the predictions of the lysis- 
only model given by Eq. ([SJ. 
Let 





(9a) 


<p = ; 

r 


(9b) 




(9c) 


r 


(9d) 


, m 
m = a„, — ; 


(9e) 




(9f) 


o-n = — ; 


(9g) 



We can non-dimensionalize the evolution equations ([5]) . 
Omitting the primes we obtain 

dTfi 

= m (1 — m) — (pmn (1 — m) — d^'m; (10a) 

at 



(in 

— = (pmn (1 — a„m) — a„n. 



Setting 



dm 
'dt 







dn 



(10b) 

(11a) 
(lib) 



we obtain three fixed points. The first is a trivial fixed 
point 

m = (12a) 
n = 0, (12b) 

which is stable when 

dm > 1. (13) 

The second corresponds to the phage extinction phase 

m = 1 — dm (14a) 
n^O, (14b) 

which is stable when 



< d„ < 1 - — 
an 



or 



1 < dm < 1 



< 



(1 - cfm) [1 - a„ (1 - d„ 



(15) 

(16) 
(17) 



where p is a root of 



an(j)p - 4>p + dn = 0. 



(18a) 
(18b) 

(19) 



The coexistence phase comes into existence and will be 
stable when 



dr. 
dn 



> Aan 
<l-p. 



(20) 
(21) 



The stability of the fixed points are governed by the Ja- 
cobian 



/ (1 - 2m) (1 - (t>n) - dn 
\ (pn (1 — 2a„TO) 

to the equations 



~(pm (1 — m) 
(j)m (1 - a„m) - rf„ 



(22) 



m (1 — m) — (pmn (1 — m) — dmm =^ (23a) 
(pmn (1 — a„m) — rf„n = 0. (23b) 

Thus we obtain the three-dimensional phase diagram 
plotted in Fig. [TJ The basin of attraction for the trivial 
case is not plotted. Region II is the basin of attraction for 
coexistence fixed point only while region III is that for the 
phage extinction. Region I will either go to coexistence 
or phage extinction, depending on the initial conditions. 




FIG. 1: (Color online) Three dimensional phase diagram for 
the lysis-only model. Region I depends on the initial con- 
ditions to flow to the phage extinction or coexistence fixed 
point. Region II and III are basins of attraction for coexis- 
tence and phage extinction fixed points, respectively. 
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C. Discussion 

As we can see, the bottom plane in Fig. [1] corresponds 
to the phase diagram in Weitz and Dushoff's model, 
where a„ = 1. When 

a> P (24) 

leading to 

an > 1, (25) 

there is a shift in the phase diagram with a rapid shrink- 
age of the basin of attraction for region II, where any 
initial condition flows to the coexistence phase. The 
boundary between region I and III also moves to larger 

which implies that the more the good and poor 

metabolisms differ from each other in the progeny num- 
ber, the easier the phages are driven out of the system. In 
order to see the effect of the phase shift more clearly, let 
us tune a„ — 1.3 while keeping all the other parameters 
as those in Fig. 2 (I) in Weitz and Dushoff's paper [SGj 
(Fig. [2]). When a„ = 1, there is a neutral fixed point 
for coexistence. However, such a fixed point disappears 
(Fig. |3]) when a„ — 1.3. The flow diagrams are generated 
by 4**^ order Runge-Kutta method. 



0.8 r 
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FIG. 2: (Color online) Flow diagram for a„ — 1, </!> = 5, ci„ — 
l,dm = 0.1. "x" denotes saddle points and "•" is for stable 
fixed points. 

In summary, we have obtained Weitz and Dushoff's 
model by detailed derivation from the individual or mi- 
croscopic level and found a small shift in the phase di- 
agram. Such a shift, as we see, can be observed experi- 
mentally by the onset of coexistence for the two species. 



1 r 




m 



FIG. 3: (Color online) Flow diagram for an — 1.3, (f} = 5, dn = 
l,dm ~ 0.1. "x" denotes saddle points and "•" is for stable 
fixed points. 



III. LYSOGENY-LYSIS MODEL 

A. Derivation of the population dynamics from an 
individual-level model 

Now we extend the lysis-only model above to incor- 
porate lysogeny and investigate the important role of 
lysogeny in host-phage dynamics. Now there are three 
types of organism in the community. They are "healthy" 
hosts, which have no integration of phage genes, lysogens, 
and free phages, which live outside bacteria membranes. 
We will label "healthy" hosts, lysogens and free phages 
by A, D, and B, respectively, with population sizes m, s 
and n. For the same reasons as in the lysis-only model, 
hosts and phages are thought of as being confined in dif- 
ferent layers characterized by different carrying capaci- 
ties. Hence both "healthy" hosts and lysogens are in the 
host layer with a total carrying capacity K. The empty 
sites in the host layer are denoted by E and their num- 
ber is K — m — s. In the phage layer, the empty sites are 
labelled by 4> as before. 

The incorporation of lysogens brings us more micro- 
scopic events. There are two pathways after phage infec- 
tions: lysis and lysogeny. Immediate lysis for temperate 
phages is the same as virulent ones, which has been char- 
acterized by events in the previous section. Lysogeny is 
an option only for temperate phages, which will be inves- 
tigated in detail here. First, there should be an event cor- 
responding to lysogen formation, i.e. a phage integrates 
its DNA into the genome of the host and turns itself into 
a prophage. Second, lysogens will survive, replicate and 
die as "healthy" hosts. Last, environments might trigger 
prophage induction, which lyses the lysogen and releases 
the prophages inside. In all, there are eighteen micro- 
scopic events, which are listed in Table Hill 
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TABLE III: Microscopic events in the lysogeny-lysis model. 



description symbols 



birth of host 


AE 


AAA 




DE 


A DD 


death of host due to longevity 


A 
/i 


— > rj 




n 
u 


— 7 rj 


death of host due to crowding 


A A 


— >■ Ah/ 




n n 


-i, DP 

— 7 U Hi 




AD 


h DE 




AD 


1 

h AE 


host-phage interactions 






• lysis under good metabolism 


AEB - 


% EEaB 


• lysis under poor metabolism 


AAB - 


U EAI3B 




ADB - 


U ED/3B 


■ lysogeny under good metabolism 


AEB 


A DE 


■ lysogeny under poor metabolism 


AAB 


A DA 




ADB 


A DD 


prophage induction 






• under good metabolism 


DEA 


► EEaB 


• under poor metabolism 


DDA 


* DEI3B 




DA A 


> AEPB 


death of free phage 


B 





Here 6, c, d, e, /, (?, h, k, p and q are constant reac- 
tion rates, a and /3 are phage reproduction numbers un- 
der good and poor metabolisms, respectively. Although 
prophage induction enhances the survival ability for lyso- 
gens in several ways, such as suppressing the latter's 
metabolism 18] through down-regulation [62j, for simplic- 
ity we have assumed the same birth and death rates for 
"healthy" hosts and lysogens. We have the condition 

a > /3, (26) 

as before. Furthermore, there are the following advan- 
tages under better metabolism: more successful and ef- 
fective infection (Eq. ()27ap ). greater possibility to lyse 
the host (Eq. (|27b|) ). and faster prophage release (Eq. 
(I27c|) '). Since lysis is controlled by lytic repressor CI 
dimers while lysogeny is regulated by CII, we do not ex- 
pect any special relationship between e and /, and p and 
q. These advantages can be expressed mathematically by 
the following inequalities: 

e + h> f + k, (27a) 

p>q. (27c) 

We draw events from the two layers the same way as in 
the lysis-only model and this results in the probabilities 
shown in Table |lVl 

From these events, we obtain the following evolution 
equations for all the three species after the calculations 
provided in Appendix B: 



— rm, ^1 ~ ^^"^ ~ (fiimn ^ ^ "^i) + (1 — 2ai) s] j- ; 

- dis + (j)2mn |l ^ ^ [(1 ^ '^si) m + (1 - 2a2i) s] 



dm 
~dt 



ds ( TO + s 

— — rs\\ 



dt~ \ K 



d2S |l - [(1 - 2022) TO + (1 - 022) s]| : 
^ = [(a - 1) 01 - a4>2\ mn |l - [(1 - 031) to + (1 - 2031) s]| 

o.d2S 1 1 - ;^ [(1 - 2032) m+{l- 032) s]| - d^n- (28c) 



.^i^^±iMl_il)(i^; (29a) 
d, = + (29b) 



2{e + h)iiv{l-uj) 
= KW ' 



(28a) 



(28b) 



where ^ 2p^.{l - ^ u) ^^^^^ 

^3-^" (29d) 



(29e) 
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021 
022 

031 

032 



2hfii> (1 — io) 

KW ' 
f + k 



2(e + /i)' 

A. 

2h'' 

q_ 

2{ae-hy 

J3q_ 

2ap 



(29f) 
(29g) 

(29h) 
(29i) 

(29j) 

(29k) 



We note that 



TABLE IV: Probabilities for the combinations in the 
lysogeny-lysis modeL 

combination probability 



AE 

DE 
A 
D 

AA 
DD 
AD 

AEB 
AAB 
ADB 
B 



M (1 - i.) (1 - c.) 
ti{l-y){l-oj) 

(1 - /.) (1 - ^) 



2m (K — m — s) 

K{K -1) 
2s{K -m-s) 

K{K ~1) 



m (m — 1) 



K 




s 


[s - 1) 


K 


(f -1) 




2ms 



flU (1 — Uj) 
flU (1 — Uj) 



K{K -1) 
2m {K ~ m ~ s) n 

K{K -1) W 
m (m — 1) n 

K{K-1)W 
2ms n 



K(K-1)W 



(f>2<^i; (30) 

< oi, 021, 022, 031, a32 < 1; (31) 

032 < 022- (32) 



We also notice some kind of symmetry in the correction 
terms such as "1 — oi" and "1 — 2ai". ai originates 
from the poor metabolism of hosts A, which indirectly 
downshifts the efficiency of phage infection and synthesis. 

In equation (|28ap . "ai" comes from the event AAB 

EA^B, while "2ai" is from ADB -A EDjiB. The factor 
"2" appears since "AD" is the same as "DA" . 

Considering 



a > 1, 



for example. 



a « 100 

for lambda phage [4§|, we approximate 

{a ~ - a(t>2 « a ((/)! ~ 4)2) ■ 
Hence equation (|28cp can be simplified as 



(33) 



(34) 



(35) 



a (01 - 02) j 1 - 77 [(1 - 031) TO + (1 - 2031) s] I + ad2S j 1 - [(1 - 2032) to + (1 - 032) s] 



dt 



K 



K 



d^n. (36) 



B. Results 

In this section, we explore the predictions of the 
lysogeny-lysis model given by equations (|28ap . (|28b|) and 

Let 



t' = rt; 
, a(j)iK 



(37a) 
(37b) 



, a4i2K 

02 = ; 

r 


(37c) 


4 = -; 

r 


(37d) 


d2 

dn = — ; 
r 


(37e) 


ds 

dn — 

3 


(37f) 


, TO 

TO 


(37g) 



9 



s 

n 



(37h) and omitting the primes, the equations after non- 
dimensionahzation become 

(37i) 
I 



dm 

'dt 
ds 

Itt 

dn 



= mil — m — s) — dim — (pimn [1 — (1 — ai) m — (1 — 2ai) s] ; (38a) 
= s (1 — m — s) - dis + (j)2mn [1 - (1 - a2i) m - (1 - 2a2i) s] — d2S [1 - (1 — 2022) m - (1 — 022) s] ; (38b) 
= - (j)2) mn [1 — (1 - 031) m — (1 - 2a3i) s] + d2S [1 - (1 - 2032) m - (1 - 032) s] — d^n. (38c) 



Formally, the fixed points can be solved by requiring 
that 



The third is the "healthy" host only fixed point 



dm 
'dt 
ds 
dt 
dn 
H 



(39a) 
(39b) 
(39c) 



However, we can only obtain four fixed points analyti- 
cally. The first is the trivial case for the extinction of all 
the species 



m = 0; 
s = 0; 
n = 0. 



(40a) 
(40b) 
(40c) 



m = 1 
s = 0; 
n = 0. 



di; 



The last is the lysogen extinction 



1 - 021 ' 



s = 0; 



1 — m — di 



[l-(l-ai)^ 



(42a) 
(42b) 
(42c) 



(43a) 
(43b) 
(43c) 



The second is the "healthy" host extinction fixed point 



m — 0: 



1 - di - ^2 
1 - d2 (1 - 022)' 



n = —s [f - (1 - 032) s] . 
"3 



(41a) 
(41b) 

(41c) 



whose existence requires that 



(01 - 02) (031 - 021) = rf3 (1 - a2i) ■ (44) 
The more interesting coexistence of all the three 
species is hard to solve analytically since the order of 
the equations 



m (1 — m — s) — dim — (pimn [1 — (1 — oi) m — (1 — 2ai) s] = 0; (45a) 
s (1 — m — s) — dis + (j}2mn [1 — (1 — 021) m — (1 — 2a2i) s] — d2S [1 — (1 — 2022) m — (1 — 022) s] ~ 0; (45b) 
(01 - 02) mn [1 - (1 - 031) TO - (1 - 2031) s] + d2S [1 - (1 - 2032) m - (1 - 032) s] - d^n = 0; (45c) 



is too high. Using a 4* order Runge-Kutta method, we C. Discussion 

found numerically a stable fixed point, shown in Fig. 2] 

As shown in Eq. p8|) . there are, in total, ten parame- 
ters so that the phase space is difficult to visualize. We 
have studied the general trend of the transition between 
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FIG. 4: (Color online) In the lysogeny-lysis model, a stable 
fixed point for the coexistence of all the three species. The 
parameters are (jji — 1,<^2 = 0.8, rfi = 0.5,^2 ~ 0.49, ds — 
0.1, ai — 021 — fflai ~ 0.1, 022 ~ 0,32 = 0.5. 




FIG. 5: (Color online) The population of the community with 
increasing phage mortality rate ds. ^^ms" indicates the sum 
of m and s. 



FIG. 6: (Color online) The population for the community 
with the increase in the lysis rate d2. 

the increase in the lysis rate, while the lysogen dimin- 
ishes. The peak in the phage population appears when 
there is a balance in the number of lysogens available to 
lyse and the lysis rate. When the lysis rate is beyond 
the threshold at 0.54, lysogen number falls dramatically 
and there is a proliferation of "healthy" hosts. The to- 
tal host population is roughly the same afterwards while 
the phage population upshifts a little with the increase 
in the "healthy" host available to infect but does not 
change further when the ratio between "healthy" hosts 
and lysogens converges. 

We have studied the effect of host mortality rate in 
Fig. [71 Obviously the total host population will fall 
monotonically when the hosts are more likely to die. We 
draw attention to the interesting peak in the phage pop- 
ulation. When the host mortality rate is low, the phage 
population is suppressed due to the overcrowding of the 
lysogens, which degrades the metabolism and hence the 
infection and synthesis of phages. When the host mor- 
tality rate is high, on the other hand, the phages have 
insufficient hosts to infect and their population also de- 
clines. 



phases, starting with the dependence on phage mortality 
rate ^3. In Fig. [5l it is shown that when the phage mor- 
tality rate is low, the systems flows into a "healthy" host 
extinction phase. The phage population decreases with 
increase in the phage mortality rate, which is very rea- 
sonable physically. For intermediate values of ds , there is 
coexistence for all the three species, while for large val- 
ues of ds, the only survival is "healthy" host, where all 
phages die out quickly out, leading to the extinction of 
lysogens. 

We show the trend of the population with increas- 
ing lysis rate d2 in Fig. [51 The phage prospers with 



D. Existence of a limit cycle 

We have noticed that the dynamics exhibits a limit cy- 
cle for some combination of parameters (Fig. (5]). In this 
section, we describe our numerical evidence for this as- 
sertion and present a plausible physical interpretation of 
our finding. In order to verify that it is a limit cycle in- 
stead of some unexpected slowing down near a putative 
stable or neutral fixed point, we have chosen an initial 
condition located inside the conjectured limit cycle. If 
there is, in fact, no real limit cycle, the dynamics will 
flow inwards no matter how slow it will be. However, as 
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FIG. 7: (Color online) The population for the community 
with the increase in the host mortality rate di. 




0.2 0.4 0.6 0.8 1 



S 



FIG. 9: (Color online) A limit cycle in the flow diagram 
with different initial conditions for parameters (pi = 1,02 = 
0.8, di = 0.5, d2 = 0.49, = 0.03, ai = 021 = asi = 
0.1, a22 = 032 = 0.5. The limit cycle is in a curved space. 
The blue curve initiated outside the cycle flows in while the 
red one from inside flows out. 
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FIG. 8: (Color online) A limit cycle in the flow diagram 
for different initial conditions with parameters (j>i = ^,4'2 ~ 
0.8, di = 0.5, d2 = 0.49, ds = 0.03, ai = 021 = aai = 

0.1, a22 ~ 0-32 ~ 0.5. 

we can see in Fig. |9l the flow indicated by the red curve 
flows out. Hence we have observed in the flow diagram 
an oscillation in the population for all the three species. 
If we inspect neighboring time steps, it appears that the 
convergence is slow, since the deviation from step to step 
is very small. However, on longer time scales, we can 
see that the convergence is an illusion. Moreover, tilt- 
ing the view angle, we see that the limit cycle is in some 
curved space instead of a single plane in Fig. 1101 In order 
to investigate the emergence of the limit cycle, we have 
scanned part of the parameter space. For example, there 
is a stable coexistence fixed point for c?i > 0.41 while 



2~, 



1.5- 




FIG. 10: (Color online) A limit cycle in the flow diagram 
with different initial conditions for parameters (j>i = 1,02 = 
0.8, di = 0.5, d2 = 0.49, ds = 0.03, ai = 021 = asi = 
0.1, a22 ~ 0-32 ~ 0.5. The limit cycle is in a curved space. 



01 = l,(/)2 = 0.8, d2 = 0.9,^3 = 0.048, ai = 021 = 031 = 
0.1, and 022 = 032 = 0.5. However, the above fixed point 
becomes unstable if di < 0.41 leading to the limit cycle. 
As we see it, such an oscillation for the population in 
the community is a manifestation of the role of lysogens 
(Fig. [TT]). When the population for host and phages 
are both small, the host will enjoy a boom because of 
good metabolism and little phage infection. Meanwhile 
prophages replicate with the fast reproduction of lyso- 
gens. Once the lysogeny-lysis switch is triggered, the de- 
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the parameters the parameters in our model to those 
which could characterize real experiments. In this sec- 
tion, we perform stochastic simulations using the Gille- 
spie's algorithm j63l [63| . which is a very efficient strategy 
to simulate chemical reactions. The reaction rates (b, c, 
d,e, / and g in Table HI and 6, c, d, e, /, h, k, p and q 
in Table IIIII) are interpreted as average probability rates 
for the occurrence of the corresponding reactions in line 
with the Gillespie algorithm, where the effect of draw 
probability is incorporated automatically. 

In the lysis-only model, the map between the two sets 
of parameters for reactions is 




b = bK- 


(46a) 


c = c; 


(46b) 


d = idif ; 


(46c) 


e = eK; 


(46d) 


7= IfK; 


(46e) 


9^9; 


(46f) 



FIG. 11: (Color online) Cartoon explanation for the limit cy- 
cle. When the population for host and phages are both small, 
the host will enjoy a boom because of good metabolism and 
little phage infection. Meanwhile prophages replicate with the 
fast reproduction of lysogens. Once the lysogeny-lysis switch 
is triggered, the destruction of lysogens will yield a huge virus 
burst. Then "healthy" host will encounter intensive phage in- 
fection and hence be suppressed. When most of the host die 
out, phage population shrinks quickly due to lack of infection. 
In this way, a cycle forms. 



where tilde is used to indicate the probability rates in 
the Gillespie algorithm. Since there are more degrees 
of freedom in choosing microscopic event rates, different 
stochastic simulations may map into the same mean field 
phase diagram. 

Our main interest is to explore the mean field limit 
cycle in the lysogeny-lysis model. We keep employing the 
tilde symbol to label the probability rates in the Gillespie 
sense and the map is 



struction of lysogens will yield a huge virus burst. Then 
"healthy" host will encounter intensive phage infection 
and hence be suppressed. When most of the host die out, 
phage population shrinks quickly due to lack of infection. 
In this way, a cycle forms. Integrating its DNA into the 
genome of a lysogen, a prophage is sheltered although it is 
temporarily dormant in the sense of viral infection. Such 
a stage assists prophages to survive demanding environ- 
mental conditions and provides an opportunity to resur- 
rect the population when there are abundant "healthy" 
hosts. Thus lysogens are perfec t g enetic reservoirs for 
phages for potential future burst [Tsl, |2]| . 



IV. STOCHASTIC SIMULATION 

Up to now, all the calculations above were carried 
out within the scope of mean field theory. As a next 
step, it is important to see to what extent such pre- 
dictions are disturbed by demographic fluctuations, and 
especially whether the limit cycle in the lysogeny-lysis 
model is stable. A second goal of this section is to link 



b^bK; 


(47a) 


c — c; 


(47b) 


d^ ^dK; 


(47c) 


e = eK; 


(47d) 


7= IfK; 


(47e) 


h = hK; 


(47f) 


fc = -kK: 
2 


(47g) 


p = pK; 


(47h) 




(47i) 


9 ^ 9- 


(47j) 



In Fig. 1121 we show a limit cycle observed in our 
stochastic simulations. It is broadly consistent with the 
mean field predictions, as can be noted easily by the ob- 
vious similarities between Fig. [T^] and Fig. [TH and Fig. 
1151 and Fig. (when we project the three-dimensional 
phase space onto two dimensions), whose relationship is 
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FIG. 12: A limit cycle in the phase space with parameters 
in the Gillespie algorithm b = 0.4, c = 0.1, d = 0.2, e = 1.2 x 
10"^°,/ = 1.2 X 10"^\p = 0.018,^ = 4.8 X 10"",fc = 4.8 x 
10-^\p = 0.54, and q = 0.27. 




FIG. 13: The projection of Fig. [12] onto the m-s plane. 



Eq. p7gp , P7h|) and (|37ip . As expected, we notice fluc- 
tuations in the stochastic simulation. For example, if Fig. 
[13 were shown in better resolution, we could see that the 
curve wiggled around the limit cycle. Usually fluctua- 
tion is two orders of magnitude smaller than the value it 
wiggles. Hence we conclude that the limit cycle is inher- 
ent to the model and robust to stochastic fluctuations, 
which serves to confirm the essential role of lysogens in 
stabilizing the cycling in the populations. 



FIG. 14: A limit cycle projected onto m-s plane in the mean 
field theory with parameters (pi — 1,(1)2 = 0.8, di = 0.5, d2 = 
0.9, di = 0.03, ai = a2i = 031 = 0.1, and 022 = az2 ~ 0.5. 
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FIG. 15: The projection of Fig. [12] onto the s-n plane. 



V. PARAMETERS IN THE MODEL 

Up to this point, all the parameters above or their val- 
ues we have explored are difficult to relate to experiment. 
The purpose of this section is to bridge the gap. 

The birth rate of the host h is medium-dependent. Usu- 
ally the expression of Lac proteins is highly suppressed 
by Lac repressors in a lacose-free medium to optimize 
energy investment and metabolism of the bacteria. In 
the above two models, we have categorized the death of 
the hosts to longevity and crowding. In fact, it is hard 
to mark a watershed clearly. Instead, what is observed 
is a population-dependent growth rate, which is a com- 
bined effect of h, c and d. Herein, the rate d for the 
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FIG. 16: A limit cycle projected onto the s-n plane in the 
mean field theory with parameters (j>i = ^,4'2 ~ 0.8, di — 
0.5, d2 = 0.9, da = 0.03, ai = a2i — aai = 0.1, and a22 ~ 
032 = 0.5. 



death of the host due to crowding is introduced artifi- 
cially to account for the actual population dependence. 
Thus, we are justified in assuming that the death rate 
of the host due to longevity c, which incorporates other 
physical and non-density-dependent factors, is fixed with 
the variation in host population. The growth rate for E. 
coli may drop to 0.2 h~^ at 37° C when glycolate serves as 
the carbon source but usually is in the range from 0.5 h~^ 
to 2.0 h~^ [23,[6^. The growth rate is species- and strain- 
specific, which for Pseudoalteromonas sp. strain SKA18 
(accessible number AF188330 in GenBank) [STj . for ex- 
ample, is an order of magnitude smaller. Similarly, lysis 
rate /, lysogeny rate k, prophage induction rate q, and 
replicate number per capita /3, which are all under poor 
metabolism, are introduced manually to characterize the 
population-dependent feature of the interactions in or- 
der to leave their population-independent counterparts 
e, h, p and a fixed. In the case of virulent phages, such 
as one in the family Siphoviri daemq attacking Pseudoal- 
teromonas sp. strain SKA18[57|, corresponding to the 
lysis-only model, the reported lysis rate spans from 0.2 
to 2.0 h~^ subject to the growth rate of the bacteria so 
that we can estimate e to be on the order of 1.0 h~^ and 
/ to be an order of magnitude smaller than e. 

For temperate phages in the lysogeny-lysis model, the 
spontaneous lysis rate is far smaller, being of the order of 
10~^ to 10"'' per generation per cell[67|. The percentage 
of lysogens is assayed through prophage induction by the 
addition of mitomycin C, UV radiation or other environ- 
mental conditions that may inhibit lambda phage repres- 
sors. Under good metabolism the lysogeny rate h for A 
phage infecting E. coli and prophage induction rate are 
on the order of 1 h""'^ and 2 h~^, respectively [68}. Their 
counterparts under poor metabolism are estimated to be 
one or two orders of magnitude smaller. For instance, the 
prophage induction rate for log-phase marine lysogens [l^ 



is on the order of 0.03 h . Replicate number per capita 
a is about 100 for phage X^^, and may be up to 600 for 
phage W-14[73|, while /3 is about 20 or 30 for both. The 
death of free phage is quite rare, which may result from 
the cleavage by proteins and depends on physical condi- 
tions such as temperature, humidity and pH values. C. 
D. Jepson and J. B. March f7ll| reported that phage A is 
highly stable, whose half life in suspension ranges from 
2.3 days at 4.2°C to 36 days at 20°C. Even if we take 
the half life be one day, the corresponding death rate g 
is on the order of 10~^ per second and can be suppressed 
by cooling down. Actually, the loss of free phage in na- 
ture, to a great extent, is through diffusion since bacteria 
are more immobile due to their large particle size com- 
pared to that of phages. In laboratory, the death rate can 
be manipulated through continuous dilution and washing 
out, and a wide range of death rates can be realized. 

When all the parameters are tuned properly, the limit 
cycle in the lysogeny-lysis model is observable in exper- 
iment. We estimate the period of the limit cycle to be 
on the order of days. Take Fig. [T3] as an example. A 
cycle there is composed of about 10, 000 computational 
steps, in other words 10, 000 events, which corresponds 
to about 120 [time unit] in the simulation. In Fig. [T3] 
the birth rate is 0.4 [time unit]^^, while in the real world 
the life cycle of an E. coli in good laboratory conditions, 
for example, is about half an hour, which is 2 hour 
Hence the cycle is 120 x 0.4/2 — 24 hours, which is one 
day 



VI. CONCLUSION 

We have derived the mean field population dynamics 
for host-phage communities both without and with lyso- 
gens. In the lysis-only model, we successfully obtained 
a description similar to the starting point assumed by 
Weitz and Dushoff , and we found that the phase dia- 
gram was modified only slightly to the difference in good 
and poor metabolism. In the lysogeny-lysis model, we 
identified the asymptotic states, which included not only 
coexistence and extinction fixed points, but population 
cycling of all microbes, lysogens and phagess. Our find- 
ings support the notion that lysogens act as a reservoir 
and are in principle amenable to experimental verifica- 
tion. We simulated the stochastic process using the Gille- 
spie algorithm and verified the robustness of our results 
to fluctuations, and especially demonstrated the stability 
of the limit cycle. 

Although complicated, our model inevitably makes 
some drastic assumptions, in addition to the most se- 
vere of all — the omission of spatial structure. In par- 
ticular, we treat "healthy" hosts and lysogens in the 
same way regarding their natural birth, death and crowd- 
ing effect. However, experimentally, the expression of 
prophage genes and the control of host gene expression by 
viral genes seem to impart to lysogens econominization in 
their metabolism[18|. When unnecessary metabolic ac- 
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tivities are suppressed, lysogens optimize their energy ex- 
penses and therefore gain some survival advantage com- 
pared to "healthy" hosts in unfavorable conditions, which 
suggests that the natural birth, death and crowding ef- 
fects of lysogens are distinct from those of "healthy" 
hosts. Hence our model is a minimal model that can 
capture the non-trivial role of lysogens in the population 
dynamics of microbe-phage communities, in addition to 
the usual predator-prey interactions. 

This work can be extended in several ways, but per- 
haps the most interesting are those which relate to the 
evolution of the field of genes distributed amongst the 
microbes, viruses and lysogens. Lysogens are genome 
carriers of not only microbes but also prophages, capable 
of yielding virus bursts when triggered by environmen- 
tal stress. In this way, the role of lysogens and viruses 
as a reservoir of genes is mediated through phage infec- 
tion and the lysogeny-lysis switch by the metabolism of 
the host. The metabolism of the host is, in turn, to 



a great extent influenced by environmental conditions. 
Thus, this model is a starting point for ecology-mediated 
evolution. It is also useful to stress that each individual 
microbe or virus constitutes a part of another organism's 
environment. Thus, the effects which our work begins 
to treat, represent a microcosm of the intricate interplay 
between ecology and evolution in microbe-virus commu- 
nities. 
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Appendix A: Transition matrices for the lysis-only 
model 



Here we provide the transition matrices, which are the 
probabilities for the change in the population in each time 
step in the lysis-only model. 



2to (K — m) ~ 

T {m + I, n\m, n) — b^il — v) (1 — oj) — = bra 

K [K — 1) 



(Al) 



26^(1 - !^) (1 -w) 26/1 (1 - J/) (1 - w) 



K - 1 



K 



(A2) 



rr^ I -,1 \ /-, . m , , ,m(m — 1) _ -r f m 1 

T (m - 1, n\m, n) = c/x (1 - v) uj— + dn{l- v) (1 - uj) _ ^ cm + dm I — - — 



cm 



— TYl 

d—- (A3) 



c/i (1 — i^) w 
K ' 



(A4) 



j_ dix{l-v){l-u:) d^i{l-v){l-u) 



2m(K^m) n ^ f m\ 

T (m — l,n + a — 1 m, n) = env (1 — cj) ^ = emn 1 ; (Ao) 

^ ' I ' y A* V ' K{K-l) W \ kJ ' ^ ' 

~_ 2eiJ.v{l-uj) 2eiiv{\-uj) 

^ {K-l)W ~ KW ■ ^ ' 
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m (m — 1) n -rn^n 
T{m-l,n + P-l\m,n)^ f^iv{l-uj) /-^; 

7_ /mJ^ (1 - _ //xz^ (1 - w) 



{K-l)W 



KW 



(A8) 



(A9) 



T {171,71- l|m, n) = 5 (1 - z^— = gn; 



(AlO) 



(All) 



I 



All the other transition matrixes are zero. Noting that that the time evolution for the probability with to hosts 
all the events in Table Hare Markov processes, we know and n phages at time t will be 



dt 



P (m, n,t) — T (m, n\m — l,n) P {771 — 1, n, i) + T (to, n|m + 1, n) P (m + 1, n, t) 



+ T (m — 1, n + /3 — 1|to, 71) + T (to, n — 1|to, n)]P (m — 1, n, t) . 

I 

Applying summations according to Eq. ([3]) , we will get 

d (m) 



dt 



{T (m + 1, n|TO, n)) — (T (to — 1, n|TO, n)) 

— (T (to — 1, n + a — 1|to, n)) 

— (T (to — 1, n + /? — 1|to, 7i)) 
(to) ' 



(S + .-)(„.)(i-if 

— e {m) (n 



c + d) (m) 



1-1 



I] H 

e K 



(Al3a) 



d{7i) 
dt 



{a — 1) (T (to — 1, n + a — 1|to, n)) 

+ (;3 - 1) (T (to - 1, n + /3 - l|m, n)) 
— (T (to, n — 1|to, n)) 

(/3-1)A(to)' 



Let 



(a - 1) e(TO) (n) 1 - I 1 
r = h + d\ 



{a-l)e K 

(A13b) 



(Al4a) 





- 1)P(to 


+ l,n 


+ /3-l,t) 


1, 7i\m, n) 


+ T(to- 


l,n + 


a — 1|to, 71) 

(A12) 




= e; 




(A14b) 
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= a — 1; 




(A14c) 




= c + d; 




(A14d) 


dn 


= 5; 




(Al4e) 




e 




(Al4f) 




(a 


-1)/ 
-l)e 


(A14g) 



which is Eq. we can arrive at Eq. ([S]). 



Appendix B: Transition matrices for the 
lysogeny-lysis model 



In this appendix, we provide details for the derivations 
of the lysogeny-lysis model. 

According to Tabic IIIII and Table IIVI we can obtain 
the following non-zero transition matrixes: 
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T (m + 1, s, n|m, s, n) = 6/x (1 - i/) (1 - w) ^ = &m f 1 - — ^ J ; (Bl) 

~ 2hix{\-v){\-i^) ^ 26/x(l-^)(l-^) 

6= K • ^^^^ 

T (m, s + 1, n|m, n) = &m (1 " y) (1 " c.) ^'jf^j^^^")'^ =hs[\- ^) ; (B3) 
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Ignoring fluctuations and correlations, we derive the evolution for population size is 
populations dynamics at the mean field level. The time 

I 



d{m) 
dt 
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dt 



d{n) 
dt 
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-{e + h) (to) {n)\l-^ 



'+h I 



+ h 



(5> 



(T (to, s + 1, n|TO, s,n)) — (T (to, s — 1, n|TO, s,n)) + {T {m — l,s + l,n — l|m, s, n)) 
— (T (to, s — 1, n + a|m, s, n)) — (T (to, s — 1, n + /3|to, s, n)) 



(to) + (s) 



- /i (to) (n) <^ 1 



K 



K 



{c + d){s) 



l-2.|j(m)+(l-||(.) 



(a — 1) (T (to — 1, s, n + a — 1|to, s, n)) + (/3 — 1) (T (to — 1, s, n + /3 — 1|to, s, n)) 



(B21a) 



(B21b) 



20 
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e + /i 

^ which is Eq. (|29p , and omit angle-brackets for simphcity, 

021 = ^; (B22h) Eq. (|B2T|) can be written as Eq. ([SS]) . 



